##########################################################################################
## https://sunduanchen.github.io/Scissor/vignettes/Scissor_Tutorial.html
library(data.table)
library(optparse)
library(ggplot2)
library(dplyr)
library(patchwork)
library(Seurat)
library(ggrepel)
library(parallel)

##########################################################################################

option_list <- list(
    make_option(c("--single_cell_file"), type = "character"),
    make_option(c("--cds_file"), type = "character"),
    make_option(c("--pathway_path"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    gene <- "MUC6"

    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic"
    single_cell_file <- paste0(work_dir,"/images/singleCell_MUC6/Scissor_STAD_MUC6_mutation.IM.CellRate.all.RData")
    out_path <- paste0("~/20220915_gastric_multiple/dna_combinePublic/images/singleCell_MUC6/")
    cds_file <- "~/ref/PCAWG_Elements/web_hg19/gc19_pc.cds.use.bed"
    pathway_path <- "~/ref/Pathway/"

}

##########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

out_path <- opt$out_path
gene <- opt$gene
single_cell_file <- opt$single_cell_file
pathway_path <- opt$pathway_path
cds_file <- opt$cds_file

dir.create(out_path , recursive = T)

##########################################################################################

sc_datasets <- load(single_cell_file, verbose = F)
coding_gene <- data.frame(fread(cds_file))
coding_gene <- unique(sapply(strsplit(coding_gene$V4 , "::") , "[" , 3))
## sc_dataset
##Idents函数定义要取的类别
##########################################################################################
Idents(sc_dataset) <- "celltype" 
## 取pit
sc_dataset <- subset(sc_dataset , idents=c("Pit"))

## 比较scissor + 和其它细胞的差异
##   Scissor_select[infos4$Scissor_pos] <- 1
##   Scissor_select[infos4$Scissor_neg] <- 2
sc_dataset$scissor <- ifelse( sc_dataset$scissor == 1 , "Positive" , "Other" )
sc_dataset$scissor <- factor( sc_dataset$scissor , levels = c("Positive" , "Other") , order = T )

Idents(sc_dataset) <- "scissor" 

## 不限制foldchange
## 显示所有基因
pit_deg <- FindMarkers(sc_dataset,
        ident.1 = "Positive" , ident.2 = "Other" , only.pos=F, min.pct=0, logfc.threshold=0
        )

## 提取coding基因
pit_deg_coding <- pit_deg[rownames(pit_deg) %in% coding_gene,]
pit_deg_coding$p_val_adj <- p.adjust( pit_deg_coding$p_val , method = "fdr" )

## 提取显著差异的基因
#pit_deg <- subset(pit_deg,(as.numeric(as.vector(pit_deg$p_val_adj))<0.05))
use_dat <- pit_deg_coding

##########################################################################################
#对原数据进行处理
use_dat$log10_q_Value <- -log10(use_dat$p_val_adj)
use_dat$log_FC <- use_dat$avg_log2FC
use_dat$gene <- rownames(use_dat)

#设置阈值
q_t <- 0.05
foldchange_t <- 1.5
logFC_cutoff <- log2(foldchange_t)
log10_P_Value_cutoff <- -log10(q_t)

data <- use_dat[,c('gene','log_FC','log10_q_Value')]
colnames(data) <- c('gene','log_FC','-log10_P_Value')
image_name <- paste0( out_path , "/DiffGene.Pit.tsv" )     
write.table( data , image_name , row.names = F , sep = "\t" , quote = F )

## Mut中显著高表达的基因
image_name <- paste0( out_path , "/DiffGene.Pit.high.tsv" )     
write.table( subset( use_dat , log_FC >= logFC_cutoff & log10_q_Value > log10_P_Value_cutoff ) , image_name , row.names = F , sep = "\t" , quote = F )

## 展示想展示的基因
data$gene <- ifelse( data$gene %in% c("GKN1" , "GKN2" , "DSC2" , "DSG2") , data$gene , "" )

col <- c('#376B6D','#9F353A' , '#F7C242','#BDC0BA' )
names(col) <- c('DOWN', 'UP', 'TARGET', "GREY" )

alpha <- 1

options(ggrepel.max.overlaps = 200)
plot1 <- ggplot(data = data,aes(x = log_FC,y = `-log10_P_Value`))+
        geom_point(data = subset(data,abs(log_FC)<logFC_cutoff),
                aes(size = abs(log_FC)),col = col["GREY"],alpha = alpha)+
        geom_point(data = subset(data,abs(`-log10_P_Value`)<log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff),
                aes(size = abs(log_FC)),col = col["GREY"],alpha = alpha)+
        geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC>logFC_cutoff),
                aes(size = abs(log_FC)),col = col["UP"],alpha = alpha)+
        geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC< -logFC_cutoff),
                aes(size = abs(log_FC)),col = col["DOWN"],alpha = alpha)+
        geom_point(data = subset(data,gene!=""),
                aes(size = abs(log_FC)),col = col["TARGET"],alpha = alpha)+
        theme_bw()+
        labs(x='log fold change\n(Scissor+ cells/All other cells)',y='-log10(adjusted p-value)')+
        geom_vline(xintercept = c(-logFC_cutoff,logFC_cutoff),lty = 3,col = 'black',lwd = 0.4)+
        geom_hline(yintercept = log10_P_Value_cutoff,lty = 3,col = 'black',lwd = 0.4) +
        geom_text_repel(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff), 
            fontface="bold", nudge_x = 0.65, nudge_y = 1, 
            aes(label = gene),size = 5,col = 'black' , face = 'bold') +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.x = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 10,color="black",face='bold'),
                axis.title.x = element_text(size = 14,color="black",face='bold'),
                axis.title.y = element_text(size = 14,color="black",face='bold'),
                strip.text.x = element_text(size = 12,color="black",face='bold'),
                axis.ticks.length = unit(0.2, "cm") ,
                axis.line = element_line(size = 0.5))

image_name <- paste0( out_path , "/DiffGene.Pit.pdf" )     
ggsave( image_name , plot1 , width = 6 )

## 通路富集在下一步做
if(1!=1){
    ##########################################################################################
    ## 通路富集
    fisher_result <- function(P0.05_all,bg_gene,gene_pathway,pathway_name){
        ta <- length(which(P0.05_all%in%gene_pathway))
        tc <- length(which(bg_gene%in%gene_pathway))
        tb <- length(P0.05_all) - ta
        td <- length(bg_gene) - tc
        data_fisher <- matrix(c(ta,tb,tc,td),nrow=2)
        fisher_res <- fisher.test(data_fisher)
        
        gene_P_in_pathway <- P0.05_all[which(P0.05_all%in%gene_pathway)]

        if(length(gene_P_in_pathway)==0){
            gene_paste <- ""
        }else if(length(gene_P_in_pathway)==1){
            gene_paste <- gene_P_in_pathway
        }else{
            gene_paste <- gene_P_in_pathway[1]
            for(i in 2:length(gene_P_in_pathway)){
                gene_paste <- paste0(gene_paste,",",gene_P_in_pathway[i])
            }
        }

        res <- data.frame(pathway=pathway_name,enrich_top=paste(ta,length(P0.05_all),sep="|"),enrich_all=paste(tc,length(bg_gene),sep="|"),bggene_in_pathway=tc,pathway_gene=length(gene_pathway),P=fisher_res$p,OR=fisher_res$estimate,gene=gene_paste)
        res
    }

    computPathway <- function(pathway_path = pathway_path , P0.05_all = P0.05_all , bg_gene = bg_gene , direction = direction){
        
        for( pathway_file in grep( "symbols.txt" , list.files(pathway_path) , value = T) ){

            file <- paste0(pathway_path , "/" , pathway_file)
            dat_pathway <- fread(file)

            pathway <- unique(dat_pathway$pathway)
            pathway_res <- data.frame(rbindlist(mclapply(pathway,function(x){
                # print(x)
                gene_pathway <- subset(dat_pathway,pathway==x)$gene
                #gene_pathway_procod <- gene_pathway[which(gene_pathway%in%bg_gene)]
            
                res_enrichment <- fisher_result(P0.05_all,bg_gene,gene_pathway,x)
                return(res_enrichment)
            },mc.cores=20)))
            results <- subset(pathway_res,OR>1)
            resultso <- results[order(results$P),]
            resultso$P_fdr <- p.adjust(resultso$P,method="fdr",n=nrow(resultso))

            out_file <- gsub("[.]symbols[.]txt" , "" , pathway_file)
            write.csv(resultso,paste0(out_path,"/",out_file,".",direction,".csv"),row.names=FALSE)
        }
    }

    ##########################################################################################
    ## 分高表达和低表达，分布做通路富集
    use_dat$padj <- use_dat$p_val_adj
    use_dat$log2FoldChange <- use_dat$avg_log2FC
    use_dat$Hugo_Symbol <- use_dat$gene

    highExpression_gene <- subset(use_dat , padj < q_t & log2FoldChange > log2(foldchange_t) )$Hugo_Symbol
    lowExpression_gene <- subset(use_dat , padj < q_t & log2FoldChange < log2(1/foldchange_t) )$Hugo_Symbol
    bg_gene <- use_dat$Hugo_Symbol

    ## 上调基因
    P0.05_all <- highExpression_gene
    direction <- "up"
    computPathway(pathway_path = pathway_path , P0.05_all = P0.05_all , bg_gene = bg_gene , direction = direction)

    ## 下调基因
    P0.05_all <- lowExpression_gene
    direction <- "low"
    computPathway(pathway_path = pathway_path , P0.05_all = P0.05_all , bg_gene = bg_gene , direction = direction)

    ## 可视化
    ## GOBP
    pathway_file <- "c5.all.v7.5.1"
    up_file <- paste0(out_path,"/",pathway_file,".up.csv")
    down_file <- paste0(out_path,"/",pathway_file,".low.csv")
    dat_up <- data.frame(fread(up_file))
    dat_down <- data.frame(fread(down_file))
    dat_up$Type <- "up"
    dat_down$Type <- "down"
    dat_plot <- rbind( dat_up , dat_down )
    dat_plot <- subset( dat_plot , P < 0.05 )
    dat_plot <- dat_plot[grep( "GOBP" , dat_plot$pathway ),]
    dat_plot$log10_P_Value <- -log10(dat_plot$P)

    mytheme <-  theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='left',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                legend.text = element_text(size = 8,color="black",face='bold'),
                axis.text.x = element_text(size = 10,color="black",face='bold' , angle = 45, hjust = 1),
                axis.text.y = element_text(size = 12,color="black",face='bold'),
                axis.title.x = element_text(size = 10,color="black",face='bold'),
                axis.title.y = element_text(size = 10,color="black",face='bold'),
                strip.text.x = element_text(size = 15,color="black",face='bold'),
                axis.line = element_line(size = 0.5))

    dat_plot <- subset( dat_plot , Type == "up" )

    dat_plot$pathway <- gsub( "GOBP_" , "" , dat_plot$pathway )
    dat_plot$pathway <- ifelse( dat_plot$pathway == "EPITHELIAL_STRUCTURE_MAINTENANCE" , "epithelial structure maintenace" , dat_plot$pathway )
    dat_plot$pathway <- ifelse( dat_plot$pathway == "MAINTENANCE_OF_GASTROINTESTINAL_EPITHELIUM" , "Maintenace of\ngastrointestinal epithelium" , dat_plot$pathway )
    dat_plot$pathway <- ifelse( dat_plot$pathway == "HOMOPHILIC_CELL_ADHESION_VIA_PLASMA_MEMBRANE_ADHESION_MOLECULES" , "Homophilic cell adhesion via \nplasma membrane adhesion molecules" , dat_plot$pathway )
    dat_plot$pathway <- ifelse( dat_plot$pathway == "DIGESTIVE_SYSTEM_PROCESS" , "Digestive system process" , dat_plot$pathway )
    dat_plot$pathway <- ifelse( dat_plot$pathway == "EXTERNAL_ENCAPSULATING_STRUCTURE_ORGANIZATION" , "External encapsulating\nstructre organization" , dat_plot$pathway )


    p1 <- ggplot(data = dat_plot,aes(y = OR,x = pathway,fill = log10_P_Value)) +
            scale_fill_distiller(palette = "Blues",direction = 1) +
            #facet_grid(Type~. , scales="free_y") +
            geom_bar(stat = "identity",width = 0.8) +
            theme_bw() +
            labs(y = "Odds Ratio",
            x = "pathway",
            title = "GOBP enrichment barplot") +
            mytheme

    image_name <- paste0( out_path , "/GOBP.Pit.pdf" )     
    ggsave( image_name , p1 , width = 6 )
}